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Abstract: We present a simple yet general and efficient approach to representation 
of computational meshes. Meshes are represented as sets of mesh entities of differ- 
ent topological dimensions and their incidence relations. We discuss a straightforward 
and efficient storage scheme for such mesh representations and efficient algorithms for 
computation of arbitrary incidence relations from a given initial and minimal set of inci- 
dence relations. The general representation may harbor a wide range of computational 
meshes, and may also be specialized to provide simple user interfaces for particular 
meshes, including simplicial meshes in one, two and three space dimensions where the 
mesh entities correspond to vertices, edges, faces and cells. It is elaborated on how 
the proposed concepts and data structures may be used for assembly of variational 
forms in parallel over distributed finite element meshes. Benchmarks are presented to 
demonstrate efficiency in terms of CPU time and memory usage. 
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Introduction 



^he computational mesh is a central component of any 
[software framework for the (mesh-based) solution of partial 
'differential equations. To reduce run-time and enable the 
Solution of large problems, it is therefore important that 
the computational mesh may be represented efficiently, 
both in terms of the speed of operations on the mesh or 
access of mesh data, and in terms of the memory usage for 
storing any given mesh in memory. 

It is furthermore important that the data structure for 
the representation of the mesh is general enough to har- 
bor a wide range of computational meshes. This generality 
must also be reflected in the programming interface to the 
mesh representation, to allow the implementation of gen- 
eral algorithms on the computational mesh. Many algo- 
rithms, such as the assembly of a linear system from a finite 
element variational problem may be implemented similarly 
for simplicial, quadrilateral and hexahedral meshes if the 
programming interface to the mesh representation does not 
enforce a specific interface limited to a specific mesh type. 
For example, if the entities on the boundary of a mesh (the 



facets) may be accessed in a similar way independently of 
the mesh dimension and not as edges in two space dimen- 
sions and faces in three space dimensions, one may use the 
same code to apply boundary conditions in 2D and 3D. 
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In (|13f ) , a very general and flexible representation of com- 
putational meshes is presented. The mesh is represented 
as a sieve, which is in general a directed acyclic graph with 
the mesh entities as points and directed edges describing- 
how the mesh entities are connected. In this paper, we 
take a slightly less gene ral approach but build on some of 
the concepts from (|13l). In particular, we will represent 
the mesh as a set of mesh entities (corresponding to the 
points of the sieve) and their incidence relations. We also 
acknowledge the works (@; 0) , where similar concepts are 
defined and where the importance of mesh iterators for 
expressing generic algorithms on computational meshes is 
advocated. 

The data structures and algorithms discussed in this 
paper have been implemented as a CH — h library and is 
distributed as part of DOLFIN, see ([H). DOLFIN is a 
problem-solving environment for ordinary and partial dif- 
ferential equations and is developed as part of the FEniCS 
project for the automation of computational mathematical 
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modeling, see (0;H). Interfaces to DOLFIN are available 
in the form of a C++ and a Python class library. 



connectivity in some detail and also introduce the concept 
mesh function. 



1.1 Design Goals 

When designing the mesh library, we had the following de- 
sign goals in mind for the mesh representation and its in- 
terface. The mesh representation should be simple, mean- 
ing that the data is represented in terms of basic C++ 
arrays unsigned int* and double*; it should be generic, 
meaning that it should not be specialized to say simpli- 
cial meshes in one, two and three space dimensions; and it 
should be efficient, meaning that operations on the mesh or 
access of mesh data should be fast and the storage should 
require minimal memory usage for any given mesh. Fur- 
thermore, the programming interface to the mesh represen- 
tation should be intuitive, meaning that suitable abstrac- 
tions (classes) should be available, including specialized 
interfaces for specific types of meshes as well as generic in- 
terfaces that enable dimension-independent programming; 
and it should be efficient, meaning that the overhead of 
the object-oriented interface should be minimized. 



1.2 Outline 

In the following section, we present the basic concepts that 
define the mesh representation and its interface. We then 
discuss the data structures of the C++ implementation of 
the mesh representation in DOLFIN, followed by a discus- 
sion of the algorithms used in DOLFIN to compute any 
given incidence relation from a given minimal set of inci- 
dence relations. Next, we demonstrate the programming 
interface to the mesh library. This is followed by a discus- 
sion of distributed (parallel) mesh data structures. Finally, 
we present a series of benchmarks to demonstrate the effi- 
ciency of the mesh representation and its implementation 
followed by some concluding remarks. 



2 Concepts 

The mesh representation is based on the following ba- 
sic concepts: mesh, mesh topology, mesh geometry, mesh 
entity and mesh connectivity. Each of these concepts is 
mapped directly to the corresponding component (class) 
of the implementation. 

A mesh is defined by its topology and its geometry. The 
mesh topology defines how the mesh is composed of its 
parts (the mesh entities) and the mesh geometry describes 
how the mesh is embedded in some metric space, typi- 
cally M n for n = 1, 2, 3. A mesh topology (Figure HJ may 
be specified as a set of mesh entities (the vertices, edges 
etc.) and their connectivity (incidence relations). Differ- 
ent embeddings (geometries) may be imposed on any given 
mesh topology to create different meshes, e.g., when mov- 
ing the vertices of a mesh in an ALE computation. Below, 
we discuss the two basic concepts mesh entity and mesh 




Figure 1: A mesh topology is a set of mesh entities (ver- 
tices, edges, etc.) and their connectivity (incidence rela- 
tions), that is, which entities are connected (incident) to 
which entities. 



2.1 Mesh Entities 

A mesh entity is a pair (d,i), where d is the topological 
dimension of the mesh entity and where i is a unique in- 
dex for the mesh entity within its topological dimension, 
ranging from to iV<j — 1 with Nd the number of entities 
of topological dimension d. We let D denote the maximal 
topological dimension over the mesh entities and set the 
topological dimension of the mesh equal to D. This is il- 
lustrated in Figure^ where each mesh entity is labeled by 
its topological dimension and index (d,i). 
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Figure 2: Each mesh entity of a mesh is identified with a 
pair (d, i), where d is the topological dimension of the mesh 
entity and where i is a unique index for the mesh entity 
within its topological dimension, ranging from to Nd — 1 
with Nd the number of entities of topological dimension d. 



For convenience, we also name common entities of low 
topological dimension or codimension. We refer to entities 
of topological dimension as vertices, entities of dimen- 
sion 1 as edges, entities of dimension 2 as faces, entities of 



codimension 1 (dimension D — 1) as facets and entities of 
codimension (dimension D) as cells. Thus, for a trian- 
gular mesh, the edges are also facets and the faces are also 
cells, and for a tetrahedral mesh, the faces are also facets. 
This is summarized in Table [T] 



Entity 


Dimension 


Codimension 


Vertex 





D 


Edge 


1 


D - 1 


Face 


2 


D-2 


Facet 


D- 1 


1 


Cell 


D 






Table 1: Named entities of low topological dimension or 
codimension. 



2.2 Mesh Connectivity 

We refer to the set of incidence relations on a set of mesh 
entities as the connectivity of the mesh. For a mesh of 
topological dimension D, there are (D+l) 2 different classes 
of incidence relations (connectivities) to consider. Each 
such class is denoted here by d — > d! for < d, d! < D. For 
any given mesh entity (d,i), its connectivity (d — > d')i is 
given by the set of incident mesh entities of dimension d' . 

Thus, for a triangular mesh (of topological dimension 
D = 2), there are nine different incidence relations of in- 
terest between the entities of the mesh. These are in turn 
— > (the vertices incident to each vertex), — > 1 (the 
edges incident to each vertex), . . . , D — > D (the cells inci- 
dent to each cell). 

For d > d' , the definition of incidence is evident. Mesh 
entity (d',i') is incident to mesh entity (d, i) if (d',i') is 
contained in {d, i). Thus, the three vertices of a triangular 
cell form the set of incident vertices and the three edges 
form the set of incident edges. For d < d' , we define mesh 
entity (d',i') as incident to mesh entity (d,i) if (d,i) is 
incident to (d',i'). It thus remains to define incidence for 
d = d! . For d, d! > 0, we say that mesh entity (d',i') 
is incident to mesh entity (d, i) if both are incident to a 
common vertex, that is, a mesh entity of dimension zero, 
while for d = 6! = 0, we say that (d', i') is incident to (d, i) 
if both are incident to a common cell, that is, a mesh entity 
of dimension D. 

Together, the set of mesh entities and the connectivity 
(incidence relations) define the topology of the mesh. Note 
that the complete set of incidence relations d — > d! for 
< d,d' < D may be determined from the single class 
of incidence relations D — »■ 0, that is, the vertices of each 
cell in the mesh. We return to this below when we present 
an algorithm for computing any given class of incidence 
relations from the minimal set of incidence relations D — > 
0. 



2.3 Mesh Functions 

We define a mesh function as a discrete function that takes 
a value on the set of mesh entities of a given fixed di- 
mension < d < D. Mesh functions are simple objects 
but very useful. A real-valued mesh function may for 
example be used to describe material parameters on the 
cells of a mesh. A boolean-valued mesh function may be 
used to set markers on cells or edges for adaptive refine- 
ment. Integer- valued mesh functions may be used to ex- 
press inter-connectivity between two separate meshes. A 
typical use is when a boundary mesh is extracted from a 
given mesh (by identifying the set of facets that are inci- 
dent to exactly one cell). One may then use a mesh func- 
tion to describe the mapping from the cells in the extracted 
boundary mesh (which has topological dimension D — 1) 
to the corresponding facets in the original mesh (which 
has topological dimension D). Note that mesh functions 
are discrete and are not meant to represent for example a 
pieccwisc polynomial finite clement function on the mesh. 



3 Data Structures 



The mesh representation as described in the previous sec- 
tion has been implemented as a small C++ class library 
and is available freely as part of the DOLFIN C++ finite 
element library, version 0.6.3 or higher. Each of the basic 
concepts mesh, mesh topology, mesh geometry, mesh en- 
tity, mesh connectivity and mesh function is realized by the 
corresponding class Mesh, MeshTopology, MeshGeometry, 
MeshEntity, MeshConnectivity and MeshFunction. All 
basic data structures are stored as static arrays of un- 
signed integers (unsigned int*) or floating point values 
(double*), which minimizes the cost of storing the mesh 
data and allows for quick access of mesh data. We discuss 
each of these classes/data structures in detail below. 

3.1 The Class Mesh 

The class Mesh stores a MeshTopology and a 
MeshGeometry that together define the mesh. The 
MeshTopology and MeshGeometry are independent of each 
other and of the Mesh. Although it is possible to work 
with the MeshTopology and MeshGeometry separately, 
they are most conveniently accessed through a Mesh class 
that holds a pair of a matching topology and geometry. 

3.2 The Class MeshTopology 

The class MeshTopology stores the topology of a mesh 
as a set of mesh entities and connectivities. For each 
pair of topological dimensions (d, d'),0 < d,d' < D, the 
class MeshTopology stores a MeshConnectivity represent- 
ing the set of incidence relations d — » d! . The mesh entities 
themselves need not be stored explicitly; they are stored 
implicitly for each topological dimension d as the set of 
pairs (d,i) for < i < Nd, where N c i is the number of 
mesh entities of topological dimension d. Thus, for each 



topological dimension, the class MeshTopology stores an 
(unsigned) integer Nd, from which the set of mesh entities 
{(d, 0), (d, 1), . . . , (d, Nd — 1)} may be generated. 

3.3 The Class MeshGeometry 

The class MeshGeometry stores the geometry of a mesh. 
Currently, only the simplest possible representation has 
been implemented, where only the coordinates of each ver- 
tex are stored. These coordinates are stored in a contigu- 
ous array coordinates of size nNo, where n is the geo- 
metric dimension and N is the number of vertices. 

3.4 The Class MeshEntity 

The class MeshEntity provides a view of a given mesh en- 
tity (d, i). The mesh entities themselves are not stored, but 
a MeshEntity may be generated from a given pair (d,i). 
The class MeshEntity provides a convenient interface for 
accessing mesh data, in particular in combination with the 
concept of mesh iterators, as will be discussed in more 
detail below. Thus, one may for any given MeshEntity 
access its topological dimension d, its index i and its set of 
incidence relations (connected mesh entities) of any given 
topological dimension d'. Specialized interfaces are pro- 
vided for the named mesh entities of Table [T] in the form 
of the following sub classes of MeshEntity: Vertex, Edge, 
Face, Facet and Cell. 

3.5 The Class MeshConnectivity 

The class MeshConnectivity stores the set of incidence 
relations d — >• d' for a fixed pair of topological dimen- 
sions (d,d'). The set of incidence relations is stored as 
a contiguous unsigned int array indices of entity in- 
dices for dimension d' entities, together with an auxiliary 
unsigned int array offsets that specifics the offset into 
the first array for each entity of dimension d0 The size 
of the first array indices is equal to the total number of 
incident entities of dimension d! and the size of the second 
array offsets is equal to the total number of entities of 
dimension d plus one. 

As an example, consider the storage of the set of inci- 
dence relations 2 — > 0, that is the vertices of each cell, for 
the triangular mesh in Figure G3 The mesh has two en- 
tities of dimension d = 2 and four entities of dimension 
d' = 0. Furthermore, each entity of dimension d = 2 is 
incident to three entities of dimension d! = 0. The array 
entities is then given by [0, 1, 3, 1, 2, 3] and the 
array of f sets is given by [0, 3, 6]. 



lr The storage is similar to the standard compressed row storage 
(CRS) format for sparse matrices, except that only the column in- 
dices need to be stored, not the values. Also note that the two 
arrays indices and offsets are private data structures of the class 
MeshConnectivity. The user is presented with a more intuitive in- 
terface, as will be demonstrated below. 
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Figure 3: The mesh connectivity 2 — > (the vertices of 
each cell) for this triangular mesh with two cells and four 
vertices is stored as two arrays indices = [0, 1, 3, 1, 
2, 3] and of f sets = [0, 3, 6]. 



3.6 The Class MeshFunction 

The class MeshFunction stores a single array of Nd 
values on the mesh entities of a given fixed dimen- 
sion d, and is templated over the value type. Typi- 
cal uses include MeshFunction<double> for material pa- 
rameters that take a constant value on each cell of a 
mesh, MeshFunction<bool> for cell markers that indicate 
cells that should be refined, and MeshFunction<unsigned 
int> to store inter- mesh connectivity or sub domain mark- 
ers. 

3.7 Minimal Storage 

The mesh data structures described above are summa- 
rized in Table [2j We note that the classes Mesh and 
MeshTopology function as "aggregate classes" that collect 
mesh data stored elsewhere, and that no data is stored 
in the class MeshEntity. All data is thus stored in the 
class MeshConnectivity (in the two arrays indices and 
offsets) and in the class MeshGeometry (in the array 
coordinates). Note that one MeshConnectivity object 
is stored for each pair of topological dimensions (d, d') for 
which the mesh connectivity has been initialized. 

As an illustration, consider the storage of a tetrahedral 
mesh with No vertices and A3 cells (tetrahedra) embed- 
ded in R 3 where we only store the set of incidence rela- 
tions D — >• 0. Each cell has four vertices, so the class 
MeshConnectivity stores 4A3 + A3 + 1 <~ 5A3 values of 
type unsigned int. Furthermore, the class MeshGeometry 
stores 3Ao values of type double. Thus, if an unsigned 
int is four bytes and a double is eight bytes, then the to- 
tal size of the mesh is 2OA3 + 24Ao bytes. For a standard 
uniform tetrahedral mesh of the unit square, the number 
of cells is approximately six times the number of vertices, 
so the total size of the mesh is 

(20A 3 + 24A ) b = (20A 3 + 24A 3 /6) b = 24A 3 b. (1) 

Thus, a mesh with 1,000,000 cells may be stored in just 



Data structure 


Principal data 


Mesh 


MeshTopology topology 
MeshGeometry geometry 


MeshTopology 


MeshConnectivity** 

connectivities 


MeshGeometry 


double* coordinates 


MeshEntity 




MeshConnectivity 


unsigned int* indices 
unsigned int* offsets 



Principal data 


Size 


MeshTopology topology 
MeshGeometry geometry 




MeshConnectivity** connectivities 




double* coordinates 


nN a 


unsigned int* indices 
unsigned int* offsets 


0{N d ) 
N d + l 



Table 2: Summary of mesh data structures. 



24 Mb. Note that if additional mesh connectivity is com- 
puted, like the edges or facets of the tctrahcdra, more mem- 
ory will be required to store the mesh. 



4 Algorithms 



In this section, we present the algorithms used by the 
DOLFIN mesh library to compute the mesh connectivity 
d — > d! for any given < d,d' < D. We assume that we 
are given an initial set of incidence relations D — > 0, that 
is, we know the vertices of each cell in the mesh. 

The key to computing the mesh connectivities of a mesh 
is to compute the connectivities in a particular order. For 
example, if the vertices are known for each edge in the 
mesh (1 — > 0), then it is straightforward to compute the 
edges incident to each vertex (0 — > 1) as will be explained 
below. The computation is based on three algorithms that 
are used successively in a particular order to compute the 
desired connectivity. As a consequence, the computation 
of a certain connectivity d — >• d! might require the com- 
putation of one or more other connectivities. We describe 
these algorithms in detail below. An overview is given in 
Figure 2] 

4.1 Build 

Algorithm [T] (Build) computes the connectivities D — > d 
and d ->• from D -» and D — > D for < d < D. In 
other words, given the vertices and incident cells of each 
cell in the mesh, Algorithm [T] computes the entities of di- 
mension d of each cell and for each such entity the vertices 
of that entity. Thus, if d = 1, then the edges of each cell 
and the vertices of each edge are computed. 

The notation of Algorithm [1] requires some explanation. 
As before, we let (d — > d')i denote the set of entities of 
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Figure 4: The three basic algorithms for computing con- 
nectivity. From the top: Build (computing connectivity 
D — > d and d — > from D — >• and £)—>£)), Transpose 
(computing connectivity d — > d! from dl — > d) and Inter- 
section (computing connectivity d — > d' from d — > d" and 
d"^d'). 



dimension d! incident to entity (d, i): 

(d ->• d')i = {(d',j) : (d'j) incident to (d,i)}. (2) 
Algorithm [T] also uses the operation 



local(_D,z) 

a > 0, 



(3) 



which denotes the set of vertex sets incident to the mesh 
entities of topological dimension d of a given cell (D, i). To 
make this concrete, consider a triangular mesh (for which 

D = 2) and take d = 1. If V 



local(Z?,i) 

a ■— > 0, then 



Vi denotes the set of vertex sets incident to the edges of 
triangle number i. The set Vi consists of three sets of 
vertices (one for each edge) and each set Vi € Vi contains 
two vertices. In addition, Algorithm [1] uses the operation 



mdex((D,j),d,Vi), 



(4) 



which denotes the index of the entity of dimension d in the 
cell (D,j) which is incident to the vertices Vi. 

We may now summarize Algorithm [1] as follows. For 
each cell (D,i), we create a set of candidate entities of di- 
mension d, represented by their incident vertices in the set 
Vi. This operation is local on each cell and must be per- 
formed differently for each different type of mesh. We then 
iterate over each cell incident to the cell (D,i) and check 
for each candidate entity S Vi if it has already been 
created by any of the previously visited cells, making sure 
that two incident cells agree on the index of any common 
incident entity. 

Algorithm 1 Build(d), computing D — > d and d —> from 

D -± and D -> D for < d < D 

k = 

for each (D, i) 

V = d ^—4 

for each {D,j) e (D — > D)i such that j < i 

local(Dj) 
Vj = d > 

for each v, G Vi 

if Vi e Vj 

I = mdcx((£>, j), d, Vi) 
{D^d) t = {D^d)iU{d,l) 

else 

(D -> d)i = {D-> d)i U (d, k) 
(d -> 0) fc = v t 
k = k + l 



4.2 Transpose 

Algorithm [2] (Transpose) computes the connectivity d — » 
d' from the connectivity dl d for d < d'. For each 
entity of dimension d' , we iterate over the incident entities 
of dimension d and add the entities of dimension d' as 
incident entities to the entities of dimension d. We may 
thus compute for example the incident cells of each vertex 
(the cells to which the vertex belongs) by iterating over the 
cells of the mesh and for each cell over its incident vertices. 



Algorithm 2 Transpose^, d'), computing d — > d' from 

d' -» d for d < d' 

for each (d',j) 

for each (d, i) € (d' 

(d -> d')i = {d- 



■>d)j 

d')iU(d',j) 



4.3 Intersection 

Algorithm [3] (Intersection) computes the connectivity d — > 
d! from d — >• d" and d" — > d' for d > d'. For each entity 
(d, z) of dimension d, we iterate over each incident entity 
(d", k) of dimension d" and for each such entity we iterate 
over each incident entity (d', j) of dimension d'. We then 
check if either (d,i) and (d',j) are entities of the same 
topological dimension or if (d',j) is completely contained 
in [d, i) by checking that each vertex incident to (d',j) is 
also incident to (d,i), in which case (d',j') is added as an 
incident entity of entity (d, i) . 

Here, d" must be chosen according to the definition of 
incidence given above. For example, we may take d" = 
to compute the connectivity D — > D (the incident cells of 
each cell) by iterating over the vertices of each cell and for 
each such vertex iterate over the incident cells. 

Algorithm 3 Intersection(d, d', d"), computing d — > d' 

from d -» ri" and d" -> d ; for d > d! 

for each (d, i) 

for each (d",fc) e (d ->• d") 4 

for each (d',j) £ (d" ->■ d') k 
if (d = d' and i ^ j) or 

(d > d' and (d' -> 0)j C (d -> 0),) 



(d — »• = (d^d') 4 U(d',j) 



4.4 Successive Application of Build, Transpose 
and Intersection 

Any given connectivity d — > d' for < d, d! < D may be 
computed by a successive application of Algorithms HHS] in 
a suitable order. In Algorithm2J we present the basic logic 
for a successive and recursive application of the three basic 
algorithms Build, Transpose and Intersection to compute 
any given connectivity. 

We illustrate this in Figure [5] for computation of the 
connectivity 2 — > 2, the incident faces of each face, for a 
tetrahedral mesh. From the given connectivity D —> 0, we 
first compute the connectivity — > D by an application 
of Transpose. This allows us to compute D D by an 
application of Intersection. The connectivity 2 — > (and 
D — > 2) may then be computed by an application of Build. 
We then apply Transpose to compute — s- 2 and finally 
Intersection to compute 2—^2. 
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Figure 5: Computing connectivity 2 — > 2 (the faces inci- 
dent to any given face) by successive application of Trans- 
pose, Intersection, Build, Transpose and Intersection. 



Algorithm 4 Connectivity^, d'), computing d — > d! by 
application of Algorithms [1] [3] 

if N d = 

Build(d) 
if N d , = 

Build(d') 
if d -> d! ^ 

return 

if d < d' 

Connectivity {d' , d) 
Transpose^, d') 

else 

if d = and d! = 

d" = D 

else 

d" = 
Connectivity(cZ, d") 
Connectivity (d", d') 
Intersection^, d', d") 



4.5 Memory Handling 

For each of Algorithms [1] O memory usage may be con- 
served by running each algorithm twice; first one round 
to count the number of incident entities, which allows the 
static data structures discussed above to be preallocatcd, 
and then another round to set the values of the incident 
entities. Furthermore, memory usage may be conserved by 
clearing incidence relations that get computed as byprod- 



ucts of Algorithms Q]-[3] when they are no longer needed. 



5 Interfaces 



In this section, we briefly describe the user interface of the 
DOLFIN mesh library. We only describe the C++ inter- 
face, but note that an (almost) identical Python interface 
is also available. 

5.1 Creating a Mesh 

A mesh may be created in one of three ways, as illustrated 
in Figure |H1 Either, the mesh is defined by a data file in 
the DOLFIN XML formaid, or the mesh is defined vertex 
by vertex and cell by cell using the DOLFIN mesh edi- 
tor, or the mesh is defined as one of the DOLFIN built-in 
meshes. Currently provided built-in meshes include trian- 
gular meshes of the unit square and tetrahedral meshes of 
the unit cube. 



// Read mesh from file 
Mesh meshOC'mesh.xml") ; 

// Build mesh using the mesh editor 
Mesh meshl; 
MeshEditor editor; 

editor . open (meshl , "triangle", 2, 2); 
editor. initVertices(4) ; 
editor. addVertex(0, 0.0, 0.0); 
editor. addVertex(l, 1.0, 0.0); 
editor. addVertex (2, 1.0, 1.0); 
editor. addVertex (3, 0.0, 1.0); 
editor . initCells (2) ; 
editor. addCelKO, 0, 1, 2); 
editor. addCelKl, 0, 2, 3); 
editor . close () ; 

// Create simple mesh of the unit cube 
UnitCube mesh2(16, 16, 16); 



Figure 6: A DOLFIN mesh may be defined either by an 
XML data file, or explicitly using the DOLFIN mesh ed- 
itor, or as a built-in predefined mesh. The last two ar- 
guments in the call to MeshEditor: :open() specify the 
topological and geometric dimensions of the mesh respec- 
tively. 



5.2 Mesh Iterators 

Mesh data may be accessed directly from the mesh, but 
is most conveniently accessed through the mesh iterator 
interface. Algorithms operating on a mesh (including Al- 
gorithms [TH3]) may often be expressed in terms of itera- 
tors. Mesh iterators can be used to iterate either over the 
global set of mesh entities of a given topological dimen- 
sion, or over the locally incident entities of any given mesh 

2 A conversion script dolf in-convert is provided for conversion 
from other popular mesh formats (including Gmsh and Mcdit) to 
DOLFIN XML format. 



entity. Two alternative interfaces are provided; the gen- 
eral interface MeshEntitylterator for iteration over en- 
tities of some given topological dimension d, and the spe- 
cialized mesh iterators Vertexlterator, Edgelterator, 
Facelterator, Facetlterator and Celllterator for it- 
eration over named entities. Iteration over mesh entities 
may be nested at arbitrary depth and the connectivity (in- 
cidence relations) required for any given iteration is auto- 
matically computed (at the first occurrence) by the algo- 
rithms presented in the previous section. 

A MeshEntitylterator (it) may be dereferenced 
(*it) to create a MeshEntity, and any member func- 
tion MeshEntity: :foo() may be accessed by it->foo(). 
A MeshEntitylterator may thus be thought of as a 
pointer to a MeshEntity. Similarly, the named mesh 
entity iterators may be dereferenced to create the cor- 
responding named mesh entities. Thus, dereferencing a 
Vertexlterator gives a Vertex which provides an in- 
terface to access vertex data. For example, if it is a 
Vertexlterator, then it->point() returns the coordi- 
nates of the vertex. 

The use of mesh iterators is demonstrated in Figure [7] for 
iteration over all cells in the mesh and for each cell all its 
vertices as illustrated in Figure [8] For each cell and each 
vertex, we print its mesh entity index. We also demon- 
strate the use of named mesh entity iterators to print the 
coordinates of each vertex. 



// Iteration over all vertices of all cells 

unsigned int D = mesh. topology O .dim() ; 

for (MeshEntitylterator c(mesh, D) ; Ic.endO; ++c) 

{ 

cout « "cell index = " << cell->index() << endl; 
for (MeshEntitylterator v(*c, 0); Iv.endO; ++v) 
{ 

cout << "vertex index = " << v->index() << endl; 

// Iteration over all vertices of all cells 
for (Celllterator c(mesh); Ic.endO; ++c) 
{ 

cout « "cell index = " << c->index() << endl; 

for (Vertexlterator v(*c) ; Iv.endO; ++v) 

{ 

cout << "vertex index = " << v->index() << endl; 

cout << "vertex coordinates = " << v->point() << endl; 

} 

} 



Figure 7: Iteration over all vertices of all cells in a mesh, 
using the general iterator interface MeshEntitylterator 
and the specialized iterators Celllterator and 
Vertexlterator. 



5.3 Direct Access to Mesh Data 

In addition to the iterator interface, all mesh data may 
be accessed directly. Thus, one may obtain an array of 
the vertices of all cells in the mesh directly from the mesh 




the mesh is ordered. Meshes may be ordered by a call to 
Mesh : : order () . 



topology, and one may obtain the vertex coordinates of the 
mesh directly from the mesh geometry. This illustrated in 
Figure H] where the same iteration as in Figure [7J is per- 
formed without mesh iterators. 



MeshTopologyfe topology = mesh. topology () ; 
MeshGeometryfe geometry = mesh. geometry () ; 
unsigned int D = topology . dim() ; 

MeshConnectivity& connectivity = topology(D, 0); 

for (unsigned int c = 0; c < topology . size (D) ; ++c) 
{ 

cout << "cell index = " << c << endl; 

unsigned int* vertices = connectivity (c) ; 

for (unsigned int i = 0; i < connectivity . size(c) ; ++i) 

{ 

unsigned int vertex = vertices [i] ; 

cout << "vertex index = " << vertex << endl; 
cout << "vertex coordinates = " 

<< geometry. point (vertex) << endl; 

} 

} 



Figure 9: Iteration over all vertices of all cells in a mesh and 
direct access of mesh data corresponding to the iteration 
of Figure [7J and Figure |U 



5.4 Mesh Algorithms 

In addition to the computation of mesh connectivity as 
discussed previously, the DOLFIN mesh library provides a 
number of other useful mesh algorithms, including bound- 
ary extraction, uniform mesh refinement, adaptive mesh 
refinement (in preparation), mesh smoothing, and reorder- 
ing of mesh entities. 

Figure [10] demonstrates uniform refinement and bound- 
ary mesh extraction. When extracting a boundary mesh, 



it may be desirable to also generate a mapping from the 
entities of the boundary mesh to the corresponding entities 
of the original mesh. This is the case for example when as- 
sembling the contribution from boundary integrals during 
assembly of a linear system arising from a finite element 
variational formulation of a PDE. One then needs to map 
each cell of the boundary mesh to the corresponding facet 
of the original mesh. (Note that the cells of the bound- 
ary mesh are facets of the original mesh.) In Figure ITTl 
we demonstrate how to extract a boundary and generate 
the mapping from the boundary mesh to the original mesh. 
The mapping is expressed as two MeshFunctions, one from 
the vertices of the boundary mesh to the corresponding 
vertex indices of the original mesh and one from the cells 
of the boundary mesh to the corresponding facet indices 
of the original mesh. 



// Refine mesh uniformly twice 
mesh.ref ine() ; 
mesh.ref ine() ; 

// Extract boundary mesh 
BoundaryMesh boundary (mesh) ; 

// Refine boundary mesh uniformly 
boundary . r ef ine ( ) ; 

// Save boundary mesh to file 
File f ile("boundary .xml") ; 
file << boundary; 



Figure 10: Uniform refinement, boundary extraction 
and uniform refinement of the boundary mesh using the 
DOLFIN mesh library. Note that the extracted boundary 
mesh is itself a mesh and may thus for example be refined. 



MeshFunction<unsigned int> vertex_map; 
MeshFunction<unsigned int> cell_map; 

BoundaryMesh boundary (mesh, vertex_map, cell_map) ; 



Figure 1 1 : Extraction of a boundary mesh and generation 
of a pair of mappings from the vertices of the boundary 
mesh to the indices of the corresponding vertices of the 
original mesh and from the cells of the boundary mesh to 
the indices of the corresponding facets of the original mesh. 



6 Parallel Considerations 



We discuss here how the concepts and data structures dis- 
cussed above can be used to assemble a global sparse fi- 
nite element operator (typically matrix) over a mesh dis- 
tributed over several processors. It is demonstrated below 
that we may reuse the concepts introduced above to dis- 
tribute the mesh. In particular, each processor owns a 



separate piece of the global mesh, which can be stored as 
a regular Mesh. Furthermore, each processor knows which 
facets of the local mesh are incident with which facets on 
other processors and this information can be stored as a 
pair of MeshFunctions. 

6.1 Simple Distribution of Mesh Data 

Let a mesh T be given and assume that the mesh has been 
partitioned into n disjoint meshes {T}™^ 1 that together 
cover the computational domain Q, C R" . Such a partition 
can be computed using for example SCOTCH, see (fl7| ). or 
Metis, see (|llt ll2l). On each processor pi, i = 0, 2, . . . , n—1, 
we store its part of the global mesh as a regular mesh and 
in addition two mesh functions Si and Ti over the facets 
of Ti. 

We thus propose to store a distributed mesh T on a set 
of processors as the set of tuples {(Ti, Si, T)}^^ where 
one tuple (Ti,Si,T) is stored on each processor pi for i = 

0,1,. ..,71-1. 

The mesh function Si maps each facet / to an integer 
j = Si(f) which indicates which (other) subdomain/mesh 
Tj that the facet / is (physically) incident with, 

Si : (D - 1, [0, 1, ... , AVi - 1]) -»• [0, 1, . . . ,n - 1]. (5) 

Here, (D — l, [0, 1, ... , N l D _ 1 — 1]) indicates that the domain 
of Si is the set of tuples {(D — l,k)} where < k < 
N}j_ 1 — 1 and N}~ ) _ 1 is the number of facets of %. Thus, 
if j = Si(f) for some j ^ i, then the facet / is shared with 
the mesh Tj. If the facet / is only incident with % itself, 
then we set <%(/) = i- 

The mesh function Ti maps each facet entity / to an in- 
teger J-i(f) which indicates which facet /' = (D — l, Ti(f)) 
of Tj for j = Si(f) that the facet / is incident (identical) 
to. If Si(f) = i, then / is not shared with another mesh 
and we set T(f) = 0. We illustrate the meaning of the 
two mesh functions Si and Ti in Figure 1121 

6.2 Parallel Assembly 

The standard algorithm for computing a global sparse op- 
erator (tensor) from a finite element variational form is 
known as assembly, see (fl9t Hoi ; 14). By this algorithm, 
the global sparse operator may be computed by assem- 
bling (summing) the contributions from the local entities 
of a finite element mesh. On each cell of the mesh, one 
computes a small cell tensor (often referred to as the "ele- 
ment stiffness matrix" ) and add the entries of that tensor 
to a global sparse tensor (often referred to as the "global 
stiffness matrix"). We shall not discuss the assembly al- 
gorithm in detail here and refer instead to SEH, but 
note that to add the entries from the local cell tensor to 
the global sparse tensor, we need to compute a so-called 
local-to- global mapping on each cell. This maps the lo- 
cal degrees of freedom on a cell (numbering the rows and 
columns of the cell tensor) to global degrees of freedom 
(numbering the rows and columns of the global tensor). 




Figure 12: A global mesh T partitioned into n — 3 meshes 
{TiYi=o ■ The mesn functions Si and Ti indicate which 
facets of % are shared with other meshes. In this example, 
we have <So((l,3)) = <?o((l,4)) = 1, indicating that facets 
(1,3) and (1,4) in To are shared with 7i- Furthermore, we 
may evaluate To at these two facets to find that the facet 
(1, 3) in 7o is incident to facet (1, 4) = To((l, 3)) in T\ and 
facet (1,4) in To is incident to facet (1,5) = J-"o((l,4)) in 
Ti. 



The assembly algorithm may be trivially parallelized by 
letting each processor pi compute and insert the cell ten- 
sors on the local mesh 77 into the global tensor. When the 
global tensor is a sparse matrix, linear algebra libraries like 
PETSc, see Q may be used to store the global tensor 
in parallel. PETSc handles the communication of matrix 
data between processors and we need only make sure that 
each processor knows how to insert entries into the global 
sparse matrix according to the local-to-global mapping of 
the global mesh. We demonstrate below that on each pro- 
cessor pi, we may (with a small amount of communication 
with neighboring processors) compute the part of the local- 
to-global mapping of the global mesh T relevant to each 
local mesh % in parallel on each processor pi, which thus 
allows us to assemble the global sparse matrix in parallel. 

6.3 Mapping Degrees of Freedom in Parallel 

In Algorithm [5J we describe how the mapping of de- 
grees of freedom may be computed on a distributed mesh 
{{%, Si,Fi)}™ =0 . 

To express this algorithm in compact form, we need to 
introduce some further notation. For each cell c in a local 
finite element mesh Ti, we assume that we can compute a 
local-to- global mapping l 1 c , which maps each local degree 
of freedom on the cell c to a global degree of freedom (for a 
numbering scheme valid on the local mesh Ti). For exam- 
ple, when computing with standard piecewise linear finite 
elements on triangles, the local-to-global mapping l 1 c may 
map a local vertex number 0, 1 or 2 on c to the correspond- 
ing global number vertex number on the mesh Ti. Thus, 
the domain of i % c is here {0, 1, 2} and the range is [0, Nq — 1], 
where Nq is the number of vertices of the mesh %■ We em- 
phasize that the local-to-global mapping is not aware of 



the global mesh T of which the local mesh % is a part. In- 
stead, it is the task of Algorithm [5] to compute (in parallel) 
a local-to-global mapping valid on the global mesh T from 
a given local-to-global mapping on each local mesh 71- 

We let Mi denote the parallel local-to-global mapping 
to be computed on each part % of the global mesh. For 
ease of notation, we express Mi as a set of tuples Mi = 
{((c, where c is a cell (number), i is the local number 

of a degree of freedom on c and / is the corresponding 
global number. The mapping Mi should be thought of as a 
function that maps a cell and a local degree of freedom (c, i) 
to the corresponding global degree of freedom /. In the 
case of standard piecewise linear elements on triangles, the 
domain of Mi is [0,N l D — 1] X {0, 1,2}, where N l D is the 
number of cells of the mesh Ti, and the range of Mi is 
[0, N — 1], where No is the total number of vertices of the 
partitioned global mesh T. We note here that the mapping 
Mi may be stored as a fixed-size array. 

To compute the parallel local-to-global mapping Mi on 
each processor, we need to iterate over the entities of the 
mesh T% and renumber the degrees of freedom. To do this, 
we introduce an auxiliary temporary mapping Mi on each 
processor that maps degrees of freedom as given by the 
local-to-global mapping on each cell c (which is only 
aware of how to map degrees of freedom internally on 71) to 
degrees of freedom as given by the parallel local-to-global 
mapping Mi (which is aware of how to map degrees of 
freedom globally on the distributed mesh). In the case 
of standard piecewise linear elements on triangles, the do- 
main of Mi is [0, Nq — 1] and the range of Mi is [0, No — 1]. 
Just as for Mi, we may think of Mi as a function but 
write it as a set of tuples (pairs) in Algorithm [5] for ease 
of notation. In a C++ implementation, Mi may be stored 
in the form of an STL map (std: :map<unsigned int , 
unsigned int>). 

Finally, we let Ni denote the number of degrees of free- 
dom on Ti not shared with a mesh Tj for j < i. This 
number can be computed on each mesh 71 from the mesh 
function Si. 

In Algorithm [5j we first let each processor pi compute 
Ni (in parallel). These numbers are then communicated 
successively from Pi-i to pi to compute an offset for the 
numbering of degrees of freedom on each mesh Ti- We 
then let each processor pi number the degrees of freedom 
(in parallel) on cells which are incident with the boundary 
of Ti and are shared with a mesh Tj for j > i. After the de- 
grees of freedom on mesh boundaries have been numbered 
on each processor, those numbering schemes are commu- 
nicated successively from pi to pj for all i < j such that 
T and Tj share degrees of freedom on a common facet. 
Finally, each processor pi may number the Ni "internal 
degrees of freedom" on %. 

The key point of this algorithm is to always let Pi number 
degrees of freedom common with pj for i < j. This num- 
bering of shared degrees of freedom is then communicated 
over shared facets, from /' to / in stage 2 of Algorithm [5] 
Since it is known a priori which facets are shared between 
two meshes Ti and Tj , one may communicate the common 



numbering for all shared degrees of freedom from pi to pj 
in one batch. 

Note that it is important that the communication of 
shared degrees of freedom in stage 2 of Algorithm [5] is 
performed sequentially, starting with processor pi receiv- 
ing the common numbering from po , then P2 receiving the 
common numbering from p and/or pi etc. This guaran- 
tees that common degrees of freedom are communicated 
from pi to all pj with i < j such that T and Tj share com- 
mon degrees of freedom, even if 71 and Tj don't share a 
common facet. For this to work, we make the assumption 
that if any two meshes % and Tj share a common degree 
of freedom on the two cells c <E T and c' <E 7j , then each 
of % and Tj must share that degree of freedom with some 
other mesh T' or Tj' respectively. In Figure Q21 we illus- 
trate this assumption (for i' = j') and give an example of 
a partition for which Algorithm [3] fails to correctly number 
all shared degrees of freedom. It is a mild assumption to 
disallow such partitions (and meshes). 

Algorithm [3] is currently not implemented in DOLFIN. 
Instead, a simple (but suboptimal) strategy where the 
computational mesh is broadcast to all processors has been 
implemented. Each processor owns a copy of the entire 
mesh and knows which part of that mesh to assemble. For 
a further discussion of the current implementation of par- 
allel assembly in DOLFIN (available with DOLFIN 0.7.2), 
see (0). 




Figure 13: Two partitions of a mesh. In the partition on 
the left, two of the meshes share only a common vertex and 
may thus share a single degree of freedom at that vertex. 
The communication of a common numbering of that degree 
of freedom is propagated over the facets incident with the 
common vertex to the neighboring mesh. In the partition 
on the right, it is not possible to propagate the numbering 
over facets (since there are no shared facets) and so Algo- 
rithm [5] will fail to compute a correct numbering scheme 
for this partition. 



Algorithm 5 {Mi} = ComputeMapping({(7I, Si, Ti)}), 
computing the local-to-global mapping in parallel for a 
mesh distributed over n processors pi, i = 0, 1, . . . , n — 1. 
— Stage 0: Compute offsets 
on each processor pi 

compute Ni 
on processor po 

offseto = 
for i = 1, 2, . . . , ?i — 1 

on processor pi 

Receive (offset^, A;_i) frompj_i 
offset* = offscti_i + Ni-i 

— Stage 1: Compute mapping on shared facets 
on each processor pi 

M i = 8 
N =_0 
h = Ni 

for each facet f E% 
3=Si(J) ' 
if j > i 

Let c G 77 be the cell incident with / 
for each local degree of freedom I on c 

if (4(0) L ) e f° r somc L 
Mi = MiU((c,l),L) 

else 

Mi = Mi U {(c,l),ki) 

M = Mu(4(Z),*i) 

ki = k% + 1 

— Stage 2: Communicate mapping on shared facets 
for j = 1, 2, . . . , n — 1 

for each facet / G Tj 
i = Sj{f) ' 
if i < j 

Receive degrees of freedom on /' from pi 
Update Afj for shared degrees of freedom 

— Stage 3: Compute mapping for interior degrees of freedom 
on each processor pi 

for each cell c E T 

for each local degree of freedom I on c 
if (t l c (0, L) G Ni for somc L 
Mi = MiU((c,l),L) 

else 

M % = Mi U {(c,l),ki) 

Ni=NiU (4(0: ki) 

k{ = ki -\- 1 



7 Benchmark Results 



In this section, we present a series of benchmarks to il- 
lustrate the efficiency of the mesh representation and its 
implementation. The new mesh library (which is available 
as part of DOLFIN since version 0.6.3) is compared to the 



old DOLFIN mesh library which is a fairly efficient C++ 
implementation, but which suffers from object-oriented 
overhead; all mesh entities are there stored as arrays of 
objects, which store their data locally in each object (in- 
cluding mesh incidence relations and vertex coordinates). 
The benchmark results were obtained on a 2.66 GHz 64- 
bit Intel processor (Q6700) running Ubuntu GNU/Linux 



for DOLFIN version 0.6.2-1 and DOLFIN version 0.7.1 re- 
spectively. 

The five test cases that are examined are the following: 
(i) CPU time and (ii) memory usage for creation of a uni- 
form tetrahedral mesh of the unit cube, (iii) CPU time for 
iteration over all vertices of the mesh, (iv) CPU time for 
accessing the coordinates of all vertices of the mesh, and 
(iv) uniform refinement of the mesh. 

In summary, the speedup was in all cases a factor 10- 
1000 and memory usage was reduced by a factor of 50 
(comparing slopes of lines in Figure IT5j) . The speedup and 
decreased memory usage is the result of more efficient algo- 
rithms and data structures, where all data is stored in large 
static arrays and objects are only provided as part of the 
interface for simple access to the underlying data represen- 
tation, not to store data themselves. Another contributing 
factor is that the old DOLFIN mesh library precomputcs 
certain connectivities (including the edges and faces of each 
cell) at startup, whereas this computation is carried out 
only when requested in the new DOLFIN mesh library, ei- 
ther as part of the iterator interface or by an explicit call 
to Mesh: :init(). We also note that from Figure HH it is 
evident that DOLFIN can be further improved in terms of 
its minimal memory footprint. 
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Figure 14: Benchmarking the CPU time for creation of 
a uniform tetrahedral mesh of the unit cube for the new 
mesh library vs. the old DOLFIN mesh library. 
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Figure 15: Benchmarking the memory usage for creation 
of a uniform tetrahedral mesh of the unit cube for the new 
mesh library vs. the old DOLFIN mesh library. 



Iteration over mesh entities 
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Figure 16: Benchmarking the CPU time for iteration over 
the vertices of each cell for the new mesh library vs. the 
old DOLFIN mesh library. 
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set of C++ classes that correspond to the basic concepts of 
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freely as part of DOLFIN jll). 
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Figure 17: Benchmarking the CPU time for iteration over 
the coordinates of each vertex for the new mesh library vs. 
the old mesh library vs. direct access of coordinate arrays 
in the new DOLFIN mesh library. 
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Figure 18: Benchmarking the CPU time for uniform re- 
finement of the unit cube for the new mesh library vs. the 
old DOLFIN mesh library. 
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